Negative thermal conductivity of chains of rotors with mechanical forcing 
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We consider chains of rotors subjected to botii thermal and mechanical forcings, in a nonequi- 
librium steady-state. Unusual nonlinear profiles of temperature and velocities are observed in the 
system. In particular, the temperature is maximal in the center, which is an indication of the 
nonlocal behavior of the system. Despite this uncommon behavior, local equilibrium holds for long 
enough chains. Our numerical results also show that, when the mechanical forcing is strong enough, 
the energy current can be increased by an inverse temperature gradient. This counterintuitive result 
again reveals the complexity of nonequilibrium states. 

PACS numbers: 44.10.+i,05.60.-k,05.70.Ln 



I. INTRODUCTION 

Thermodynamic properties of non-equilibrium station- 
ary states are poorly understood. They are usually char- 
acterized by currents of conserved quantities, such as en- 
ergy, flowing through the system. When stationary states 
are close to equilibrium states, linear response theory is 
effective and explains common macroscopic phenomena 
like Fourier's law: In a system in contact with two ther- 
mostats at different temperatures, the heat flux is propor- 
tional to the temperature gradient (as long as the relative 
difference between the two temperatures is small). 

In contrast, there is no general theory to describe sys- 
tems in a stationary state far from equilibrium, and the 
corresponding macroscopic properties seem to depend on 
the specific details of the dynamics. 

In this article, we investigate numerically the energy 
transport properties of a simple one-dimensional system, 
a chain of N rotors, in a stationary state far from equi- 
librium. Many studies considered one-dimensional chains 
of oscillators subjected to a temperature gradient [l|, Q- 
Here we consider both thermal and mechanical forcings, 
obtained as follows: The leftmost rotor is attached to a 
wall and put in contact with a Langevin thermostat at 
temperature Tl, while the rightmost rotor is subjected 
to a constant external force F and put in contact with 
another Langevin thermostat at temperature Tr. 

What we observe in our numerical experiments is that 
the combined effect of these two generalized forces can 
reduce the current instead of increasing it. This counter- 
intuitive effect is observed for large mechanical forcings 
F , when Tr is increased while Tl remains fixed (see Fig- 
ure IS] below) . The mechanical forcing induces a negative 
current (from the right to the left). When Tr is increased 
while Tl remains fixed, one would naively expect that the 
(negative) thermal forcing is also larger, and thus that 
the negative current should be (in absolute value) larger. 
In contrast to this expectation, we observe that, in this 



case, the current is reduced. This strange effect does not 
appear if, instead, Tl is lowered and Tr is fixed (see Fig- 
ure |6] below) , in which case the current indeed becomes 
larger in absolute value. 

We are unable to provide explanations to the above 
described phenomena. We believe that such behaviors 
show the complexity of non-equilibrium stationary states 
far from equilibrium, and also suggest that Fourier's law 
is only valid close to equilibrium. A naive extension of 
the definition of thermal conductivity to genuinely non- 
equilibrium settings can give negative values to this quan- 
tity (whence the somehow provocative title of this arti- 
cle). 

In the sequel of this article, we first describe the system 
we work with, and the numerical integrator we have used 
(see Sec. nil). We next turn to studying various properties 
of the system. In particular, we numerically check that 
local equilibrium holds for systems large enough, despite 
the fact that, globally, the system is out of equilibrium 
(see Sec. IIIIBj) . In Sec. IIII (Jl we study how the current 
depends on the magnitude of the mechanical force and 
on the temperatures that are imposed on both ends of 
the chain. All these numerical studies are performed for 
chains of increasing lengths. 



II. DESCRIPTION OF THE SYSTEM 

The configuration of the system is described by the 
positions (angles) q — {qi, . . . ,qN) of the rotors, which 
belong to the one-dimensional torus 27rT, as well as their 
associated (angular) momenta p = {pi, . . . ,pn). The 
masses of the particles are set to 1 for simplicity. The 
Hamiltonian of the system is 



N 



Hiq,p)^Yl 



i^ + {l-cosr. 



(1) 



where we have set Vi — qi — qi-i for i > 2 and ri = qi. 

We consider a system with free boundary conditions 
on the right end, whose evolution equations read: 



dq^ ^Pidt, 



dpi = f sin(g,;+i - qi) - s\\\{q^ - qi^i) ]dt, i^ 1, N, 
dpi = (^sin(q2 ~ qi) - sm{qi)jdt 



IPi 



dt + y/¥{Kdwl, 



dpN = [F - sin(qjv - qN-i)]dt 



N 



-lPNdt + ^2-iT^dWi 

(2) 
where W^ and W^^ are independent standard Wiener 
processes, and 7 > determines the strength of the cou- 
phng to the thermostat. In the sequel, we work with 
7 = 1. Note that the external constant force F is non- 
gradient since it does not derive from a periodic potential. 
We checked the robustness of the results we describe 
below with respect to the choice of boundary conditions. 
We indeed also considered fixed boundary conditions on 
the right end (this amounts to adding an extra force 
— sin(gjv) to the last atom). In particular, we checked 
that our counterintuitive results on the behavior of the 
thermal current as a function of the strength of the non- 
gradient force are still observed with these boundary con- 
ditions. 



A. Equilibrium and nonequilibrium states 

If _F = and Tl = Tr, = T , the system is in equi- 
librium, and the unique stationary measure is given by 
the Gibbs measure at temperature T associated with the 
Hamiltonian Q. When Tl ^ Tr, with F = 0, the prop- 
erties of the non-equilibrium stationary state have been 
studied numerically by various authors |3|-|6[. The ther- 
mal conductivity of the system, defined as the station- 
ary energy current multiplied by the size of the system 
and divided by the temperature difference Q , has a finite 
limit for large system sizes, even though the rotor chain is 
a momentum conserving one-dimensional system |8|, |9[ . 
Besides, as the average temperature T increases above 
the value 0.5, the thermal conductivity decreases dra- 
matically [3|. 

If -F 7^ 0, the system is out-of-equilibrium even if 
Tl = Tr (recall indeed that F is non-gradient). The 
force, in the stationary state, induces an energy current 
towards the left. The stationary state cannot be com- 
puted explicitly and, if F is large, linear response theory 
cannot be used to obtain information about the conduc- 
tivity of the system. 

If Tl < Tr, there are two mechanisms that separately 
generate an energy current towards the left of the system: 
The mechanical force F and the thermal force given by 
the temperature gradient. It seems however difficult to 



separate the contributions of each mechanism. The nu- 
merical experiments reported below show that these two 
mechanisms arc not necessarily additive, and that one 
mechanism may reduce the effect of the other one, lead- 
ing to counterintuitive results. 



B. Numerical integration 

The numerical integration of ^ is performed using a 
splitting strategy where the Hamiltonian part of the evo- 
lution is integrated with the Verlet scheme [ifll. The 
fluctuation-dissipation parts, with the additional non- 
gradient force, are Ornstein-Uhlenbeck processes and can 
thus be integrated analytically. We have thus used the 
following algorithm: 



pT 


= ap^ + ai^G^, 


p'n 


= F + a{pl-F)+aKGl, 


p": 


= K, ^y^hN, 


P^" 


^' 2 ^q^^^ '-^ '' 


ir' 


= ,^+Atpr^/^ 


p^' 


_ „+i/2 Aiaif +1 „+i/2^ 

^' 2 dq^^ '^ 


where a = 


3xp(-7Af), fiL = V{l-a^)n 



■\/(l - a^)TR, and H is given by (P). In turn, G" and 
G^ are independent normal Gaussian random variables. 
Recall also that the friction parameter 7 is set to 1. The 
three first lines of the above algorithm consist in exactly 
integrating the Ornstein-Uhlenbeck processes on pi and 
Pn, whereas the three last lines arc based on the standard 
Verlet algorithm. 

The time-step At = 0.05 ensures that the energy con- 
servation in the Verlet scheme is accurate enough. While 
there might be some time-step bias in the value of the 
currents, the qualitative conclusions are robust with re- 
spect to the choice of the time-step. 



III. PROPERTIES OF THE NONEQUILIBRIUM 

SYSTEM 



This section is organized as follows. First, we discuss 
the existence of a stationary measure for the dynam- 
ics ^. Under the assumption that such a stationary 
measure exists, we establish some relations that are con- 
sistent with physical intuition. We next point out that 
this system shows some very surprising features. For in- 
stance, the temperature profile is non-monotonic and a 
maximum is observed in the center of the system, while 
the velocity profiles are very nonlinear. Despite these 
nonlocal features, we show that local equilibrium holds. 
We finally turn to investigating the dependence of the 
stationary energy current on F, Tl and Tr. 



A. Stationary measure 

We believe that there exists a unique smooth station- 
ary measure for the dynamics ^. However, as far as we 
know, there is no rigorous resuh in this direction for ro- 
tor chains, even in the case F = 0. Indeed, the standard 
techniques (see for instance [ll|, [l2| ) used to prove exis- 
tence and uniqueness of an invariant measure for chains 
of oscillators under thermal forcing do not apply here. 

A possible pathology for rotor chains is that the (in- 
ternal) energy concentrates locally on one or several ro- 
tors, which rotate faster and faster. Since the interaction 
forces arc bounded, it may not be possible to prevent 
this fast rotation. In practice, we have not observed 
such catastrophes in the parameter regime we consid- 
ered, but the kinetic temperature profiles presented in 
Figure [T] (obtained from the variance of the momenta, 
with the previously mentioned caveat on the interpreta- 
tion of this quantity) are quite unexpected and show that 
the internal energy tends to be larger in the middle of the 
chain. This picture also allows to understand what hap- 
pens when the imposed temperatures at the right and 
left ends change: The maximal temperature in the chain 
is almost unchanged, but the position of the maximum is 
displaced. This shows that the linear response correction 
to the stationary measure is necessarily nonlocal. Such 
nonlocal effects were already observed in nonequilibrium 
exclusion processes [ij, [Ijl . 
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i{pn) for all i > 2, while Fi :— (sinri) = F — j{{pn) + 

iPi))- 

The balance between the average work done by the 
force and the energy dissipated by the thermostats is 
given by 

= F{pn) + lin - (P?)) + 7(Tr - {P%)), (3) 

as can be seen by noticing that the average variation of 
the total energy H is zero. Moreover, the entropy pro- 
duction inequality (obtained by computing the variations 
of the relative entropy with respect to the invariant mea- 
sure, see e.g. [l^l) gives 

T{:\n - (pD) + T^\Tn. - {p%)) < 0. 

In the case Tl = Tr, = T, this relation, combined 
with (|3]), yields F{pn) > 0. Therefore the stationary 
momentum on the right end has the same sign as the 
driving force, as expected. 

Figure [5] shows that the momentum profile is not linear, 
and that its derivative is maximal where, according to 
Figure [1] temperature is maximal. We also observe on 
Figure [3] (top) that the profile seems to become steeper 
in the thermodynamic limit. 




200 



400 



600 
site 



800 



1000 



FIG. 2: Average momenta for chains of length N = 1024, 
with F = 1.6. 



FIG. 1: Kinetic temperature profiles for chains of length A'^ 
1024, with F = 1.6 and different temperature gradients. 



Some interesting relations can nonetheless be obtained 
under the assumption that the stationary state exists. 
We denote by (•) the expectation with respect to the sta- 
tionary measure. First, a constant profile of force settles 
down in the bulk. Taking expectations in ^ indeed gives 



(sinrj+i) = (sinrj), 

(sinr2) = (sinri) +7(^1), 
(sinr-Af) ^ F -j(pn). 



t^l,N, 



This leads to the following profile: Fi := (sinri) = F 



B. Local equilibrium and thermodynamic limit 

A very interesting question is whether nonequilibrium 
systems are locally close to equilibrium. This issue was 
considered in [1^ for systems subjected to thermal forc- 
ings only. We check the local equilibrium assumption in 
three steps. 

(i) We study the agreement between the local kinetic 
temperature (defined as the variance of the veloci- 
ties) and the local potential temperature. The lat- 
ter is obtained as follows. First, we numerically 



precomputc the function 



V{r)cxp{~V{r)/T)dr 
exp{-V{r)/T) dr 

which, to a given temperature, associates the canon- 
ical average of the potential energy V{r) = 1 — cos r 
of one bond. The local potential temperature at 
bond i is then defined as the value Ti such that 
g{Ti) is equal to the time average of the potential 
energy of the bond r^ along the trajectory defined 
by © . The results presented in Figure [3] (bottom) 
show that the two local temperatures are quite dif- 
ferent for small systems, but are identical for larger 
ones. Besides, as the length of the system increases, 
the profiles become more symmetric. 




temperature is maximal (since this is the location 
where the disagreement between the local kinetic 
and potential temperatures is the strongest). The 
results presented in Figure|3]show that the empirical 
distributions of p and r at the site imax are in ex- 
cellent agreement with the Gibbs distributions with 
the same parameters (average velocity Pj, tempera- 
ture Ti), namely 



and 



Z,-:iexp[-(p-p,)V(2T,)] dp 



Zp-iexp hy(r)/T,]dr, 



except for the smallest systems (say, N < 512). 
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FIG. 3: Rescaled profiles for systems of increasing size N = 2'' 
with k = 7, . . . , 14: the x variable is the site index i divided 
by A''. The value of the nongradient force is F — 1.6 and 
Tl — Tr. = 0.2. Top: momenta. Bottom: kinetic (solid lines) 
and potential (dashed lines) temperatures. 



(ii) We check that the individual distributions of p and r 
are in accordance with a local Gibbs equilibrium. 
To this end, we build the histograms of the mo- 
menta and distances at the site imax where the local 



data histogram 
empirical density 
reference density 




1 2 3 4 5 6 

distance 

FIG. 4: Top: Empirical distribution of momenta at the site 
Jmax (where the temperature is maximal), and comparison 
with the local Gibbs equilibrium with the same average and 
variance. Bottom: Empirical distribution of the distances at 
bond imax and comparison with the local Gibbs equilibrium 
with the same average energy. Both plots correspond to a 
chain of length N = 1024, with F = 1.6 and Tl = Tr = 0.2. 



(iii) We check that momenta and distances are inde- 
pendent. To this end, we compare the joint law 

V' = ^('^imaxi-Pimax) ^f (j't^.x 7 K„ax ) ^nd the product 
law obtained from the tensor product of the indi- 
vidual distributions of these two variables (denoted 
respectively by ■0r(''imax) ^^.nd ippiPi^^^))- More 
precisely, for a given number n of sample points 
(obtained by subsampling a long trajectory every 
10"' steps), we check that the distance 



[0,27r]xl 



^"(r,p)-Vr(?-)V'p(p) drdp (4) 



between these two distributions indeed decreases as 
the inverse square-root of the number of configura- 
tions used to build the histograms. Again, this is 
true for systems large enough. Figure [5] shows that 

5n -^ n-1/2. 
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FIG. 5: Decrease of the error 5„ defined by (Q as a function 
of the number of sample points ?i. Estimated rate of decrease: 
Sn ~ 28.7 X n-o-^^l 



C. Behavior of the energy current 

We consider the following situations: 

(i) same temperatures on the left and on the right: 
(Tl,Tr) = (0.20,0.20) or (0.15,0.15); 

(ii) hot left end and cold right end: (Tl,Tr) — 
(0.25,0.15), (0.20,0.15) or (0.25,0.20); 

(iii) cold left end and hot right end: (TLiTpj) — 
(0.15,0.25). 




FIG. 6: Comparison of the currents with fixed tempera- 
ture on the right end and increasing temperatures on the 
left end. From top to bottom: decreasing system sizes 
iV = 2048, 1024, 512, 256, 128 (the ordering is the same for 
all situations considered; for the longest systems, we have 
considered forces < F < 1.6; for the shortest ones, we have 
considered the range F £ [0; 2.4]). 



Currents are computed as a function of the magnitude F 
of the non-gradient forcing term for systems of different 
lengths: N = 128, 256, 512, 1024, 2048. RecaU that local 
equilibrium holds at the leading order, so that the energy 
current is induced by the first order corrections in 1/N. 
We first compare the currents when the temperature 
on the right end is fixed (see Figure |6]). As expected, 
the negative current induced by the mechanical forcing 
is reduced by the opposite, positive thermal current. 

We next compare the currents at fixed temperature dif- 
ference Tr — Tlj for different average temperatures (see 
Figure [7]). In this case, we observe that, for strong me- 
chanical forcings, the current is enhanced when the av- 
erage temperature decreases, while the opposite happens 
when the mechanical forcing is small. 

We finally turn to the most interesting situation. The 
temperature on the left end is fixed and the tempera- 
ture on the right end varies (see Figure [8]). In this case, 
counterintuitive results are observed for large mechanical 
forcings: The total current is enhanced as Tr decreases. 
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FIG. 7: Comparison of the currents with fixed temperature 
difference. From top to bottom: decreasing system sizes 
iV = 2048, f 024, 5f 2, 256, 128 (the ordering is the same for 
all situations considered). 



even though, in such a situation, the thermal gradient is 
in the opposite direction. The mechanical forcing induces 
a negative current, while the thermal gradient induces (in 
the absence of any force) a positive current. The com- 
bined effect of both mechanical and thermal forcings in- 
duces a negative current larger (in absolute value) than 
the one in the absence of any thermal gradient! 



IV. DISCUSSION OF THE RESULTS 

In conclusion, for large mechanical forcings F ^ wc ob- 
serve that 

(a) when Tr is fixed, the current varies qualitatively as 
when there is no mechanical forcing: The absolute 
value of the current increases when Tl decreases, 
which means that the current induced by the ther- 
mal forcing and the current induced by the mechan- 
ical forcing are somewhat additive. In this case, a 
positive thermal conductivity is observed (for a fixed 
value F of the mechanical forcing, considering only 
the response in the limit when Tr — Tl — > 0). 
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FIG. 8: Comparison of the currents for a fixed temperature at 
the left end and various temperature difi^erences. From top to 
bottom: decreasing system sizes A^ = 2048, 1024, 512, 256, 128 
(the ordering is the same for all situations considered). 



(b) when Tl is fixed, the current has a surprising behav- 
ior: Its absolute value increases when Tr decreases. 
This means that the thermal forcing, which is naively 
expected to reduce the current induced by the me- 
chanical forcing, actually enhances it. In this case, a 
negative thermal conductivity is observed (again, for 
a fixed value F of the mechanical forcing). 

A possible interpretation is based on the fact that, for 
such a system, the thermal conductivity is a decreasing 



function of the temperature when F is large (see Fig- 
ure [7]). It is possible that, by lowering Tr and thus in- 
creasing the conductivity at the right end, one makes the 
system more sensitive to the mechanical forcing. The 
increased mechanical current may hence counterbalance 
the increased opposite thermal current. 

An interesting question which wc did not discuss here is 
the scaling of the energy current as a function of the sys- 
tem size when F ^ 0. Some preliminary results suggest 
that the thermal conductivity is finite, as when F = 0, 
but this question definitely calls for additional studies. 
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